#########################################################################################
# Replication for: Juen et al. "Compulsory Vaccination", Politics. Forthcoming.
#########################################################################################

# LOAD PACKAGES

library(tidyverse)
library(prediction)
library(haven)
library(Rfast)
library(interflex)

# Load Data

df <- read_spss("qualtrics_data.sav")

# Treatment Coding

df$pflicht_treatment <- NA
df$pflicht_treatment[!is.na(df$pflicht_base)] <- 1
df$pflicht_treatment[!is.na(df$pflicht_foe)] <- 2
df$pflicht_treatment[!is.na(df$pflicht_friend)] <- 3
df$pflicht_treatment[!is.na(df$pflicht_rki)] <- 4
df <- filter(df, !is.na(pflicht_treatment)) # deletes incomplete responses

df$pflicht_treatment <- factor(df$pflicht_treatment,
                               levels = 1:4,
                               labels = c("No cue",
                                          "Out-party cue",
                                          "In-party cue",
                                          "Expert cue"))

# Best Party Evaluation 

df$best_eval <- apply(df[,paste0("party_eval_",1:6)],1,max)

# Second Best Party Evaluation

df$second_best_eval <- apply(df[,paste0("party_eval_",1:6)],1,function(x){
  n <- length(x)
  sort(x,partial=n-1)[n-1]
})

# Difference between First and second best party

df$best_diff <- df$best_eval - df$second_best_eval

# Outcome

df$support_pflicht <- NA
df$support_pflicht[!is.na(df$pflicht_base)] <- df$pflicht_base[!is.na(df$pflicht_base)] %>% as_factor %>% as.numeric()
df$support_pflicht[!is.na(df$pflicht_foe)] <- df$pflicht_foe[!is.na(df$pflicht_foe)] %>% as_factor %>% as.numeric()
df$support_pflicht[!is.na(df$pflicht_friend)] <- df$pflicht_friend[!is.na(df$pflicht_friend)] %>% as_factor %>% as.numeric()
df$support_pflicht[!is.na(df$pflicht_rki)] <- df$pflicht_rki[!is.na(df$pflicht_rki)] %>% as_factor %>% as.numeric()

# Populism

df$populism <- as.numeric((df$populism_1_1 + 
                             df$populism_1_2 + 
                             df$populism_1_3 + 
                             df$populism_2_1 + 
                             df$populism_2_2 + 
                             df$populism_2_3)/6)

# Pol Interest in 3 Categories

df$polint3 <- ifelse(df$pol_int < 3, "Low",
                     ifelse(df$pol_int == 3, "Moderate", "High")) %>%
  factor(levels = c("Low", "Moderate", "High"))

# Left Right

df$left_right <- as.numeric(df$left_right)

# Party Choice

df$afd_vote <- as.numeric(df$party_vote == 3)
df$left_vote <- as.numeric(df$party_vote == 6)
df$populist_vote <- as.numeric(df$afd_vote | df$left_vote)

# Party vote

df$party_vote2 <- as_factor(df$party_vote)

####################################################################################
# OLS Models 
####################################################################################

mod1 <- lm(support_pflicht ~ pflicht_treatment + populism + left_right, data = df)
mod2 <- lm(support_pflicht ~ pflicht_treatment + populism * left_right, data = df)
mod3 <- lm(support_pflicht ~ pflicht_treatment*best_eval + populism + left_right, data = df)
mod4 <- lm(support_pflicht ~ pflicht_treatment*best_diff + populism + left_right, data = df)
mod5 <- lm(support_pflicht ~ pflicht_treatment*pol_int + populism + left_right, data = df)
mod6 <- lm(support_pflicht ~ pflicht_treatment*party_vote2 + populism + left_right, data = df)

texreg::texreg(list(mod1,mod2),
               stars = c(0.001, 0.01, 0.05, 0.1),
               booktabs = TRUE)

texreg::texreg(list(mod3,mod4, mod5),
               stars = c(0.001, 0.01, 0.05, 0.1),
               booktabs = TRUE)

####################################################################################
# Plot Interaction LR * Populism (assuming linear interaction effect as done by lm)
####################################################################################

preds <- prediction(mod2, at = list("populism" = c(1, 3, 5),
                                    "left_right" = 1:11)) %>% 
  summary %>% 
  as.data.frame()

names(preds)[1:2] <- c("Populism", "LeftRight")

preds$Populism <- factor(preds$Populism,
                         levels = c(1, 3, 5),
                         labels = c("Low", "Moderate", "High"))

ggplot(preds, aes(x = LeftRight, 
                  y = Prediction,
                  fill = Populism)) +
  geom_ribbon(aes(ymin = lower,
                  ymax = upper), 
              alpha = .5,
              color = NA) +
  geom_line(aes(lty = Populism)) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  scale_fill_manual("Populism", 
                    values = c("#576e5f", "#815a82", "steelblue")) +
  ylab("Predicted Support for\nCompulsory COVID-19\nVaccination") +
  xlab("Left-Right") +
  scale_x_continuous(labels = c("Left", rep("",9), "Right"), breaks = 1:11)

ggsave(file = "left_right_populism.pdf",
       width = 10,
       height = 5)
ggsave(file="figure1.jpg", width = 20, height = 15, units = "cm")

#########################################################################
# 3D Plot using GAM for Appendix (based on new package interflex-syntax)
#########################################################################

dfx <- select(df,
              Support = support_pflicht,
              populism = populism,
              treat = pflicht_treatment,
              lr = left_right) %>%
  as.data.frame()

dfx$lr <- as.numeric(dfx$lr)

jpeg("gam.jpg", res=300,
     width = 30, height = 15, units = "cm")
interflex(Y="Support",
          D="populism", 
          X="lr", 
          data=dfx,
          Ylabel = "\n\nSupport\nf. Vaccination",
          Xlabel = "\nLeft-Right",
          Dlabel = "Populism",
          estimator = "gam",
          angle = c(-320, -170, 340, -80))
dev.off()

####################################################################################
# FIGURE 2
####################################################################################

preds <- prediction(mod3, at = list("pflicht_treatment" = unique(df$pflicht_treatment),
                                    "best_eval" = 1:11)) %>% 
  summary %>% 
  as.data.frame()

names(preds)[1:2] <- c("pflicht_treatment", "best_eval")
names(preds)

outparty <- filter(preds, pflicht_treatment %in% c("No cue", "Out-party cue"))
outparty$facet <- "Out-party cue"
inparty <- filter(preds, pflicht_treatment %in% c("No cue", "In-party cue"))
inparty$facet <- "In-party cue"
expert <- filter(preds, pflicht_treatment %in% c("No cue", "Expert cue"))
expert$facet <- "Expert cue"

newdata <- bind_rows(outparty, inparty, expert)
newdata$D <- ifelse(newdata$pflicht_treatment == newdata$facet,
                    "Treatment", "Control")

ggplot(newdata, aes(x = best_eval, 
                  y = Prediction,
                  fill = D)) +
  geom_ribbon(aes(ymin = lower, 
                  ymax = upper,
                  alpha = D), 
              color = NA) +
  geom_line(aes(lty = D)) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Predicted Support for\nCompulsory COVID-19\nVaccination") +
  xlab("Level of 'Best Party' Evaluation") +
  scale_alpha_manual("",values = c(0.3, 0.7)) +
  scale_fill_manual("",values = c("steelblue", "orange")) +
  scale_linetype_manual("",values = c(2,1)) +
  scale_x_continuous(labels = c("Very\nbad", 
                                rep("",9), 
                                "Very\ngood"), 
                     breaks = 1:11) +
  facet_wrap(~ facet)

ggsave(file = "best_party_effect.pdf",
       width = 8,
       height = 4)

ggsave(file="figure2.jpg", width = 25, height = 15, units = "cm")

####################################################################################
# FIGURE 3
####################################################################################

preds <- prediction(mod4, at = list("pflicht_treatment" = unique(df$pflicht_treatment),
                                    "best_diff" = 0:10)) %>% 
  summary %>% 
  as.data.frame()

names(preds)[1:2] <- c("pflicht_treatment", "best_diff")

outparty <- filter(preds, pflicht_treatment %in% c("No cue", "Out-party cue"))
outparty$facet <- "Out-party cue"
inparty <- filter(preds, pflicht_treatment %in% c("No cue", "In-party cue"))
inparty$facet <- "In-party cue"
expert <- filter(preds, pflicht_treatment %in% c("No cue", "Expert cue"))
expert$facet <- "Expert cue"

newdata <- bind_rows(outparty, inparty, expert)
newdata$D <- ifelse(newdata$pflicht_treatment == newdata$facet,
                    "Treatment", "Control")

ggplot(newdata, aes(x = best_diff, 
                    y = Prediction,
                    fill = D)) +
  geom_ribbon(aes(ymin = lower, 
                  ymax = upper,
                  alpha = D), 
              color = NA) +
  geom_line(aes(lty = D)) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Predicted Support for\nCompulsory COVID-19\nVaccination") +
  scale_alpha_manual("",values = c(0.3, 0.7)) +
  scale_fill_manual("",values = c("steelblue", "orange")) +
  scale_linetype_manual("",values = c(2,1)) +
  xlab("Distance of 'Best Party' to 'Second-Best Party'") +
  scale_x_continuous(labels = c("Tied\nfirst\nplace", 
                                rep("",9), 
                                "Clear\nfirst\nplace"), 
                     breaks = 0:10) +
  facet_wrap(~ facet)

ggsave(file = "best_diff_effect.pdf",
       width = 12,
       height = 5)


ggsave(file="figure3.jpg", width = 30, height = 15, units = "cm")

###########################################################
# MIP (appendix)
###########################################################

mip <- data.table::fread("mip_fgw.csv",
                         data.table = F)

mip$Date <- lubridate::dmy(mip$Date)

mip <- mip[,!grepl("V[0-9]+", names(mip))]

mip <- pivot_longer(mip,
                    cols = names(mip)[2:7],
                    names_to = "Topic",
                    values_to = "Percent")

mip$Year <- lubridate::year(mip$Date)
mip$Month <- lubridate::month(mip$Date)

mip %>%
  filter(Year == 2020) %>%
  group_by(Topic, Year, Month) %>%
  summarise(Percent = mean(Percent, na.rm = TRUE)) -> mip

mip <- mip %>%
  mutate(MIP = case_when(grepl("^Aus", Topic) ~ "Migration/Refugees",
                         grepl("^Bil", Topic) ~ "Education",
                         grepl("^Cor", Topic) ~ "COVID-19",
                         grepl("^Rent", Topic) ~ "Pensions",
                         grepl("^Soz", Topic) ~ "Social Inequality",
                         grepl("^Umw", Topic) ~ "Environment/Climate"))

ggplot(mip, aes(x = Month,
                y = Percent, color = MIP)) +
  geom_line() +
  geom_point() +
  theme_bw(base_size = 14) +
  xlab("Month of 2020") +
  ylab("'Most Important Problem'\n(in %)") +
  scale_x_continuous(breaks = 1:12,
                     label = month.abb) +
  theme(legend.position = "bottom")

ggsave("mip.pdf", width = 10, height = 5)
ggsave(file="figure_mip.jpg", width = 30, height = 15, units = "cm")

####################################################################################
# Figure 4
####################################################################################

preds <- prediction(mod5, at = list("pflicht_treatment" = unique(df$pflicht_treatment),
                                    "pol_int" = 1:5)) %>% 
  summary %>% 
  as.data.frame()

names(preds)[1:2] <- c("pflicht_treatment", "pol_int")

ggplot(preds, aes(x = pol_int, 
                  y = Prediction,
                  color = pflicht_treatment)) +
  geom_pointrange(aes(ymin = lower, ymax = upper),
                  position = position_dodge(width = 0.3)) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  guides(color = guide_legend("", ncol = 2)) +
  ylab("Predicted Support for\nCompulsory COVID-19\nVaccination") +
  xlab("Political Interest") +
  scale_x_continuous(breaks = 1:5,
                     labels = c("Low", "...", "Moderate", "...", "High"))

ggsave(file = "pol_int_interaction.pdf",
       width = 8,
       height = 5)
ggsave(file="figure4.jpg", width = 20, height = 15, units = "cm")

####################################################################################
# APPENDIX PARTY 
####################################################################################

preds <- prediction(mod6, at = list("pflicht_treatment" = unique(df$pflicht_treatment),
                                    "party_vote2" = unique(df$party_vote2))) %>% 
  summary %>% 
  as.data.frame()

names(preds)[1:2] <- c("pflicht_treatment", "party_vote2")

preds$party_vote3 <- case_when(grepl("^Bünd", preds$party_vote2) ~ "Green Party",
                               grepl("^Die", preds$party_vote2) ~ "Left Party",
                               grepl("^Eine ", preds$party_vote2) ~ "Other Party",
                               grepl("^Ich wü", preds$party_vote2) ~ "Non-Voters",
                               TRUE ~ paste(preds$party_vote2))

ggplot(preds, aes(x = pflicht_treatment, 
                  y = Prediction,
                  color = pflicht_treatment)) +
  geom_pointrange(aes(ymin = lower, ymax = upper),
                  position = position_dodge(width = 0.3)) +
  theme_bw(base_size = 14) +
  theme(legend.position = "bottom") +
  ylab("Predicted Support for Compulsory COVID-19 Vaccination") +
  xlab("") +
  facet_wrap(~ party_vote3, scales = "free") +
  coord_flip() +
  guides(color = F)

ggsave(file = "party_interaction.pdf",
       width = 10,
       height = 7)
ggsave(file="figure_A1.jpg", width = 25, height = 15, units = "cm")
